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Abstract: A critical point (E) is expected in QCD on the temperature (T) versus baryonic 
chemical potential (//) plane. Using a recently proposed lattice method for \i ^0 we study 
dynamical QCD with nf =2+1 staggered quarks of physical masses on Lt = 4 lattices. Our 
result for the critical point is Tg = 162 ± 2 MeV and \xe = 360 ± 40 MeV. For the critical 
temperature at fi = we obtained T c = 164 ± 2 MeV. This work extends our previous 
study [Z. Fodor and S.D.Katz, JHEP 0203 (2002) 014] by two means. It decreases the light 
quark masses (m Ut d) by a factor of three down to their physical values. Furthermore, in 
order to approach the thermodynamical limit we increase our largest volume by a factor of 
three. As expected, decreasing m u ^ decreased he- Note, that the continuum extrapolation 
is still missing 
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I. Introduction. 

QCD is an asymptotically free theory, thus its high temperature and high density phases 
are dominated by partons (quarks and gluons) as degrees of freedom rather than hadrons. 
These phases play an important role in the particle physics of the early universe, of neutron 
stars and of heavy ion collisions (for a clear theoretical introduction and a review see Q). 

Extensive experimental work has been done with heavy ion collisions at GSI, CERN 
and Brookhaven to study the strong interactions at high temperatures and non-vanishing 
baryon densities (the latter correspond to non-vanishing baryonic chemical potentials) and 
to explore the fi-T phase diagram. It is a long-standing open question, whether a critical 
point (E) exists on the \x-T plane, and particularly how to predict theoretically its location 

II, |H . At this point E the phase transition is of second order and long wavelength fluctua- 
tions appear, which results in characteristic experimental consequences, similar to critical 
opalescence. Passing close enough to (he,Te) one expects simultaneous appearance of sig- 
natures. The observables exhibit non-monotonic dependence on the control parameters ||, 
since one can miss the critical point on either of two sides. 

The location of this critical point is an unambiguous, non-perturbative prediction of 
the QCD Lagrangian. Unfortunately, until recently no ab initio, lattice analysis based on 
QCD was done to locate the endpoint. Models with infinitely large strange quark mass 
(m s =oo) were used (e.g. ||), suggesting that he ~ 700 MeV. The result is sensitive to 
the strange quark mass (he should be smaller for smaller m s ). For realistic cases these 
techniques can not predict the value of he even to within a factor of 2-3. 

Lattice QCD at non-vanishing baryon density should, in principle, give an unam- 
biguous answer. Though QCD at finite fi can be formulated on the lattice |J, standard 
Monte-Carlo techniques can not be used at fj, 7^ 0. The reason is that for non-vanishing 
real \i the functional measure -thus, the determinant of the Euclidean Dirac operator- 
is complex. This fact spoils any Monte-Carlo technique based on importance sampling. 
Several suggestions were studied earlier to solve the problem. Unfortunately, none of them 
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was able to give the phase line or locate (/j>e,Te)- About two years ago new techniques 
appeared, with which moderate chemical potentials could be reached on the lattice. 

In two recent papers we proposed a new method Jj], ||] to study lattice QCD at finite 
T and \x. The idea was to produce an ensemble of QCD configurations at fi=0 and at 
the corresponding transition temperature T c (or at any other physically motivated point 
for which importance sampling works). Then we determined the Boltzmann weights Q 
of these configurations at /u 7^ and at T lowered to the transition temperatures at this 
non-vanishing \i. An ensemble of configurations at a transition point was reweighted to an 
ensemble of configurations at another transition point. With this technique a much better 
overlap was observed than by reweighting pure hadronic ensemble to a transition one [10]. 
We illustrated the applicability of the method in nj=4 dynamical QCD Q and in nj =2+1 
dynamical QCD ||. The phase line with the critical end-point and the equation of state 
1 11, 12] were determined. 

A less CPU demanding, truncated version of the overlap-improving multi-parameter 
reweighting was also applied using an improved lattice action. The approach can be sum- 
marized as follows. Instead of evaluating the Boltzmann weights exactly one can expand 
it in the chemical potential and use the first terms of the expansion. This modified tech- 
nique was also able to give the phase diagram [JD]] and the equation of state [14] (estimates 
based on derivative extrapolation for the chiral endpoint was also reported in conference 
proceedings [15, [l(|). The success of the overlap-improving multi-parameter reweighting 
was analyzed in Ref. |T^| . 

A completely independent method is based on the the fact that at imaginary chemical 
potential importance sampling works. The result on the phase digram can be analytically 
continued to real chemical potentials. This technique gave the phase line both in nj=2 
I IS ] and in nj=A [IS], which are consistent with the results of the overlap-improving multi- 
parameter reweighting method [^, ^, Analytic continuation was also used to estimate 
the location of the critical endpoint p^] for nj=3. 

Recently several other new techniques were suggested (see e.g. |2l], ^]), which will 
be most probably tested in the near future. The recent developments of lattice QCD at 
non- vanishing chemical potentials are reviewed by Refs. |24][ . 

In this paper we use our original suggestion 0, |8| and evaluate the reweighting Boltz- 
mann factors exactly. We determine the volume (V) dependence of the zeros of the partition 
function on the complex gauge coupling {(3) plane. Based on this volume dependence we 
determine the type of the transition as a function of \i. The endpoint //£ is given by 
the value at which the crossover disappears and finite volume scaling predicts a first order 
phase transition. These finite T calculations are done on Lt = 4 lattices. In order to set the 
physical scale we determine the pion, kaon and rho masses (m n , mx, mp), and the Sommer 
pq] scale (Ro) at T=0. Our quark masses are realistic, the strange quark mass and the 
light quark masses are set about to their physical values. Having determined the lattice 
spacing we transform our result to physical units and give T c , the location of (^e,Te) and 
show the phase diagram separating the hadronic phase and the QGP. 

The present work is a significant improvement on our previous analysis by two means. 
We increased the physical volume by a factor of three and decreased the light quark masses 
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by a factor of three, down to their physical values. Due to these improvements the compu- 
tational effort of the present work was 140 times larger than that of the previous analysis 
||. As expected, decreasing the light quark masses resulted in a smaller he- Increasing the 
volumes did not influence the results, which indicates the reliability of the finite volume 
analysis. In order to give the final answer to {he,Te) the most important step remained 
is the continuum extrapolation. Note, however, that for the present physical problem de- 
creasing the lattice spacing by a factor of 2 increases the CPU costs by approximately three 
orders of magnitude. 

The remaining part of the paper is organized as follows. In Section 2 we summarize the 
overlap-improving multi-parameter reweighting technique and the method of the Lee- Yang 
zeros, which can be used to separate the crossover and first order transition regions. Section 
3 contains the details of the T^O and T=0 simulations. Those who are not interested in 
the lattice details could directly move to Section 4, in which we present our results and 
conclude. 

2. Overlap-improving multi-parameter reweighting at /i ^ 0. 

The partition function of lattice QCD with nj degenerate staggered quarks (for an intro- 
duction see e.g. |2(|) is given by the functional integral of the bosonic action Sf> at gauge 
coupling [3 over the link variables U, weighted by the determinant of the quark matrix M, 
which can be rewritten j7j as 



where m is the quark mass, /i is the quark chemical potential and rif is the number of 
flavours. For non-degenerate masses one uses simply the product of several quark matrix 
determinants on the 1/4-th power. Standard importance sampling works and can be used 
to collect an ensemble of configurations at m w , j3 w and [i w (with e.g. Re([i w )=0 or non- 
vanishing isospin chemical potential). It means we treat the terms in the curly bracket 
as an observable -which is measured on each independent configuration- and the rest as 
the measure. By simultaneously changing several parameters e.g. (3 and (i one can ensure 
that even the mismatched measure at (5 W and \x w samples the regions where the original 
integrand with (5 and \x is large. In practice the determinant is evaluated at some \i 
and a Ferrenberg-Swendsen reweighting |]] is performed for the gauge coupling (3. The 
fractional power in eq. fl2,l| ) can be taken by using the fact that at fi = fi w the ratio of 
the determinants is 1 and the ratio is a continuous function of the chemical potential. The 
details of the determinant calculation can be found in Ref. j8). 

In the following we keep /i real and look for the zeros of the partition function on 
the complex f3 plane. These are the Lee- Yang zeros p7|| . Their V— > oo behavior tells the 
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Table 1: T ^ and T = results. The upper part is a summary of the Lee- Yang zeros obtained 
at different chemical potentials for m Mj d=0.0092 and to s =0.25. We indicate by 6,8,10,12 and oo 
the spatial extensions -and their extrapolation- of our Lt = 4 lattices. The lower part shows the 
measured T = observables for three /3 and three m u ,d values at to s =0.25 on 12 3 • 24 lattices. 

difference between a crossover and a first order phase transition. At a first order phase 
transition the free energy oc log Z(f3) is non-analytic. Clearly, a phase transition can appear 
only in the oo limit, but not in a finite V. Nevertheless, the partition function has 
Lee- Yang zeros at finite V. These are at "unphysical" complex values of the parameters, 
in our case at complex (3-s. For a system with a first order phase transition these zeros 
approach the real axis in the oo limit (the detailed analysis suggests & 1/V scaling). 
This oo limit generates the non-analyticity of the free energy. For a system with 
crossover the free energy is analytic, thus the zeros do not approach the real axis in the 
oo limit. 

3. T ^ and T = simulations for n f =2+l. 

Using the formulation described above we study 2+1 flavour QCD at T ^ on Lt = 

4, L s = 6, 8, 10, 12 lattices with m u ^ = 0.0092 and m s = 0.25 as bare quark masses. 
Note, that these mass parameters approximately correspond to their physical values. At 
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T = we use 12 3 -24 lattices. Three different couplings (/3=5.090, 5.160, 5.190) are studied 
in order to determine the non-perturbative /^-function. Chiral extrapolation in the light 
quark masses at T=0 are done by using three different mass parameters (m nj d=0.02, 0.04, 
0.06). For each parameter set 3000 configurations were generated. For generating the field 
configurations the R algorithm is applied. The microcanonical stepsize is always set to half 
of the light quark mass. We use a modified version [28] of the MILC collaboration's code 

®. 

At T ^ we determined the complex valued Lee- Yang zeros, /?o, for different V-s as a 
function of fi. Their V— > oo limit was calculated by a (3q(V) = + £/V extrapolation. 
The results (listed in Table |) are obtained by generating 100,000; 100,000; 100,000 and 
150,000 configurations on our L s =6,8,10 and 12 lattices, respectively. The determinant 
calculation was carried out after every 50 trajectories. Thus our results are based on a few 
thousand independent configurations. 

Figure ffl shows Im(/3o°) as a function of \x enlarged around the endpoint Hend- 

The 

picture is simple and reflects the physical expectations. For small /i-s the extrapolated 
Im(/3^°) is inconsistent with a vanishing value, and the prediction is a crossover. Increas- 
ing [i the value of Im(/3^°) decreases, thus the transition becomes consistent with a first 
order phase transition. The statistical error was determined by a jackknife analysis using 
subsamples of the total L s =6,8,10 and 12 partition functions. Our primary result in lattice 
units is fj, en d = 0.1825(75). 
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Figure 1: In^/^) as a function of the chemical potential. 



Table |] contains also the T = results. To set the physical scale we used an average 
obtained from Rq (0.5 fm) and m p (770 MeV). Note, that these two quantities give some- 
what different scales (the difference is fairly small, approximately 10-15%). This difference 
is expected to dissappear when one approaches the continuum limit with physical quark 
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masses. 

Setting the scale leads to the final results of the analysis. As we already discussed, the 
quark masses, used to determine the endpoint, correspond approximately to their physical 
values. The pion to rho mass ratio, extrapolated to our parameters, is 0.188(2) (its 

physical value is 0.179), whereas the pion to K mass ratio in the same limit is 0.267(1) (its 
physical value is 0.277). 

Along the transition line (3 decreases, thus the lattice spacing increases. During the 
reweighting procedure we did not change the lattice quark masses, thus the quark masses 
changed in physical units. As a consequence, the transition line slightly deviates from the 
line of constant physics. We corrected for this effect by using our previous results at larger 
m ud fH an d by following the observation that transition lines for slightly different quark 
masses are practically parallel (see [jll], [l^] ) . 

Let us estimate the applicability of the method approaching the continuum limit. 
In the present analysis with physical quark masses the evaluation of the eigenvalues was 
somewhat less costly than the production of the configurations. Extending the analysis 
to even larger volumes, approximately upto 4-16 3 , the eigenvalue determination remains 
subdominant. Using these larger volumes reduces the error on fi en d to a level, which is 
not even needed (uncertainties due to finite lattice spacing are more important). For finer 
lattices the eigenvalue evaluation goes with and the configuration production goes at 
least with L^. According to numerical estimates |ll], [l?]] the applicability range of 
the overlap-improving multi-parameter reweighting technique, along the transition line, 
scales with [i oc V 1 with 7 ~ 1/3. (Note, that the systematic study of the lattice spacing 
dependence of 7 has not been performed yet.) In the scaling region the chemical potential 
of the endpoint is constant in physical units. It scales with V ' 25 in lattice units for fixed 
physical volumes. Taking into account the marginal difference between 7 ~ 1/3 and 0.25, 
one concludes that the determinant evaluation remains subdominant and even with the 
present technique one might successfully approach smaller lattice spacings. 

4. Results and conclusions. 

Figure ^ shows the phase diagram in physical units, thus T as a function of the baryonic 
chemical potential (which is three times larger then the quark chemical potential). 

The transition temperature at vanishing baryonic chemical potential is T c =164±2 MeV. 
Note, that this value is somewhat smaller than our previous result (T c =172±3 MeV of Ref. 
H). This is a known phenomenon: smaller quark mass results in smaller transition tem- 
perature (see e.g. |30|] ). 

The curvature of the crossover line separating the QGP and the hadronic phases is 
given by T/T c = 1 — Cfi 2 B /T 2 with C=0. 0032(1). This value is somewhat smaller than our 
previous curvature ||| or other values in the literature [H], [l4|, Note, that compared 
to other analyses, we took into account an additional effect which reduced the value of C. 

1 Note, that this difference is smaller than the difference in our previous analysis ||. This improvement 
is due to the choice of the physical light quark masses. 
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The small change of the mass parameter on the line of constant physics (caused by the 
change of the lattice spacing) slightly decreases the curvature. 

The endpoint is at Tg = 162 ±2 MeV, he = 360 ±40 MeV. As expected, he decreased 
as we decreased the light quark masses down to their physical values (at approximately 
three-times larger m u ^ the critical point was at he=720 MeV; see ||). 
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Figure 2: The phase diagram in physical units. Dotted line illustrates the crossover, solid line the 
first order phase transition. The small square shows the endpoint. The depicted errors originate 
from the reweighting procedure. Note, that an overall additional error of 1.3% comes from the 
error of the scale determination at T=0. Combining the two sources of uncertainties one obtains 
T E = 162 ± 2 MeV and he = 360 ± 40 MeV. 

The above result is a significant improvement on our previous analysis || by two 
means. We increased the physical volume by a factor of three and decreased the light 
quark masses by a factor of three. Increasing the volumes did not influence the results, 
which indicates the reliability of the finite volume analysis. Clearly, more work is needed 
to get the final values. Most importantly one has to extrapolate to the continuum limit. 
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